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I. INTRODUCTION 

A central problem in the theory of the static structure of classical liquids is to find 
a simple and efficient way to obtain the pair correlation functions from the inter-particle 
forces in pairwise interacting fluids. Exact statistical mechanics allows to write the 

formal solution of such problem as the coupled set of equations: 

1 + hij{r) = exp[-/50ij(r) + hij{r) - Cij{r) + Bij{r)] (1) 

and 

kj{r) = Cij{r) + J2p^ f dr'cu{r')hij{\r - r'|) (2) 
I 

where hij{r) and Cij{r) are the (total) and indirect correlation functions for atomic pairs of 
species i and j, pi is the number density of the l-th component and f3 = 1/kT. The functions 

n 

Bij{r), named bridge functions after their diagrammatic characterization ll| are junctionals 
of the total correlation functions, i.e. their value at distance r depends on the values of all 
the correlation functions at all distances. 

The basic difficulty with equations ((H) and (j2|) is that we do not have an explicit and 
computationally efficient relation between Bij{r) and the correlation functions, so we have to 
resort to approximations. The results of the last three decades of research have shown that 
it is possible to make progress by approximating the bridge functionals Bij{r) by functions 
of the indirect correlation functions 7jj(r) = hij{r) — Cij{r) (approximate closures). Once we 
have an explicit form for Bij{jij{r)), the resulting integral equations (H)) and ((21) , although 
approximate, can provide excellent results for the static structure of liquids. Moreover, be- 
sides the original focus on the structural properties, in recent years interest has grown toward 
using approximate integral equations to obtain thermodynamics and the phase diagrams of 
liquids and liquid mixtures [3|. 

In particular, Kjellander and Sarman |^ and Lee have derived an approximate but 
useful formula for the chemical potential of a fluid requiring only the knowledge of the 
correlation functions at the thermodynamic state of interest. Their formula is based on two 
main approximations. The first is the same assumption from which integral equations are 
derived, i.e. that the bridge functions Bij{r) are local functions of the corresponding indirect 
correlation functions. The second stronger assumption is that the only dependence of the 
bridge functions from the thermodynamic state is through the indirect correlation functions. 
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Thus, the functional dependence of Bij{'^ij) is the same for all the states. 

In this paper we want to investigate via direct numerical computer simulation the two 
approximat ions . 

Up to now, numerical studies of the bridge functions and of the accuracy of the local 

riQ 

approximation have been limited to the case of one component systems p, 13| or electrolytic 
solutions [8]. We feel that two-component systems deserve more interest for many reasons: 
i) there are strong indications that the approximate universality of the bridge functions is 
not valid in multicomponent systems, ii) the phase diagrams of multicomponent systems are 
richer and more interesting than those of pure fluids, and iii) it turns out that modeling the 
bridge functions for multicomponent systems is much more difficult than for pure systems. 

We have studied, through Monte Carlo simulation, the bridge functions of a few systems of 
non-additive hard spheres (NAHS) mixtures, including the limiting cases of additive (AHS) 
mixtures and one component system. In particular we are interested in a direct check of 
the local hypothesis for the functional relations between bridge and correlation functions in 
binary mixtures. To this aim we use the so called Duh-Haymet plots Q]. These are plots 
of the partial bridge functions Bij as a function of the partial indirect correlation functions 

lij- 

The paper is organized as follows. In section |n] we summarize the equations we used 
to evaluate the cavity correlation functions from which the bridge functions can be easily 
obtained and we provide the relevant technical details of the numerical calculations. In 
section Uni we present and discuss our numerical results. 



II. CALCULATION OF THE CAVITY AND BRIDGE FUNCTIONS 



A. Theory 

The binary NAHS system is a fluid made of hard spheres of two species. One specie, here 
named 1, with diameter Rn and number density pi and another specie (2) with diameter 
i?22 and number density with a pair interaction potential that can be written as follows 



Mr) = — , (3) 
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where R12 = {Rn + -^22)72 + a, with a being the non additivity parameter. We will also 
study various special cases as the one component system, and the binary mixture of additive 
hard spheres (AHS) a = 0. We can rewrite Eq. (JIJ to obtain the partial bridge functions 

Bab{r) = lnyab{r) - ^ab{r) , (4) 

where yab{f) are the partial cavity functions 

Vabir) = Qabir) exp[/?0ab(r)] , (5) 

here gij{r) = 1 + hij{r) are the partial radial distribution functions. Notice that both the 
cavity functions and the indirect correlation functions are everywhere continuous then also 
the bridge is so. 

In the region outside the hard cores, in a hard sphere (HS) system, the cavity correlation 
functions coincide with the pair distribution functions gij{r). In order to determine the 
relationship between the partial bridge functions and the partial indirect correlation func- 
tions within the hard cores, we need to calculate the partial cavity functions. There are two 

one which uses Henderson's equation and 
the direct simulation method of Torrie and Patey We decided to use the first method 
which is accurate at small r. 

For a binary mixture the like cavity functions can be obtained from the following canonical 
average 

VZa_ 



yaa[r l„2a, 




,ia>2 it=l 



, (6) 

Ni,N2,V,T 

where a,b = 1,2 with b a, Vi^j^ is the distance between particle i of specie a and particle 
j of specie b, = exp(/5/Xa)/A^ is the activity of specie a, Ha its chemical potential, and A 
the de Broglie thermal wavelength, V is the volume, Na the number of particles of specie 
a, so that the prefactor Vza/Na = exp^Pfi^^'^) where /i^^"^ is the excess chemical potential 
of specie a. The notation (. . .)ni,N2,v,t indicates the canonical average at fixed number of 
particles, volume and temperature. 

So to calculate Vaaij) we need to introduce in the system of Na + A*";, particles labeled 
lb, ... , Afe, 2a, ... , (A^+l)a a test particle 1^ placed a distance r from particle 2^ and calculate. 
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at each Monte Carlo step, the interaction of this particle with all the particles of the system 
except particle 2^. 

We immediately realize that when r = we must have 

^aa(O) = 1 , (7) 

since the configurations where particle 2^ overlaps with other particles of the system are 
forbidden. Moreover, by taking into account that yabi'f') = fi'afe(^) for r > Rab and from the 
asymptotic value of the partial pair distribution functions follows that 

lim yaair) = e'"'''^^^ . (8) 

r— >oo 

The unlike cavity functions can be obtained from the following canonical average 

f \ Vzi_ 



'^^'^ ' exp { —(3 



N2 Ni+l 



, (9) 

Ni,N2y,T 



J2>1 n>l 

So to calculate yi2{r) we need to introduce in the system of Ni + N2 particles labeled 
I2, • • • , N2-, 2i, . . . , (A^+l)i a test particle li placed a distance r from particle I2 and calculate, 
at each Monte Carlo step, the interaction of this particle with all the particles of the system 
except particle I2. 

Now there is no simple argument to guess the contact value of yi2. All we can say is that 
we must have 2/12(0) < 1. At large r we still have 

lim yi2(r) = e-^'^i'' . (10) 



B. Numerical implementation 

Monte Carlo simulations were performed with a standard NVT Metropolis algorithm 
using N = 4000 particles. Linked lists have been used to reduce the computational cost. 
To calculate the partial pair distribution functions we generally used 5.2 x 10^ Monte Carlo 
steps, where one step corresponds to the attempt to move a single randomly chosen particle, 
and incremented the histograms once every 20 x 4000 steps. To calculate the partial cavity 
functions we used 1.6 x 10^ Monte Carlo steps and incremented the histograms once every 
2 X 4000 steps. The acceptance ratio was adjusted to values between 10% and 40%. 
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The Monte Carlo simulation returned the gabi'f') over a range not less than 8.125i?ii for 
the densest system. In all the studied cases, the pair distribution functions attained their 
asymptotic value well inside the maximum distance they were evaluated. Thus, it has been 
possible to obtain accurate Fourier transforms of the total correlation functions [hab{k)] 
(it was necessary to cure the cusps at contact in the partial pair distribution functions by 
adding to them H{Rab — r)gab{Rab), H being the Heaviside step function, before taking 
the Fourier transform and removing its analytical Fourier transform afterwards). To obtain 
the partial indirect correlation functions we first calculated the partial direct correlation 
functions [ca{,(A;)] using the Fourier transform of the Ornstein-Zernike equation (j2I) and then 
we got the Fourier transform of the indirect correlation functions ^ij{k) = hij{k) — Cij{k) 
which is the transform of a continuous function in real space and then is safe to transform 
back numerically to obtain 7afe(r). 



III. NUMERICAL RESULTS 



We carried on simulations on the following systems: (A), one component HS; (B), equimo- 
lar binary mixture of AHS; (C), equimolar binary mixture of NAHS with equal like diameters 
and negative non additivity; (D), equimolar binary mixture of NAHS with equal like diame- 
ters and positive non additivity; and (E), equimolar binary mixture of NAHS with different 
like diameters. In all these cases we have drawn the corresponding Duh-Haymet plots, i.e. 
we plot, for each distance, the pairs (i?jj(r), 7jj(r)). 

When we are outside the hard core the partial bridge functions (jD) reduces to 

Bab{r) = lngab{r) - 7ab(r) . (11) 

and we can obtain the bridge functions directly from the pair correlation functions solving 
the OZ equation ^ to get the partial indirect correlation functions 'jab- 

To realize the Duh-Haymet plots when we are within the hard core regions, we first 
calculated the cavity functions yab as explained in section |H] and then the bridge functions 
(up to an additive constant, the excess chemical potential (3^1^'^ ) from their definition (^. 
Estimating the excess chemical potential from the long range behavior of the cavity functions 
[see equations ^ and (fTUj) ] we where able to find the full bridge functions. Since the cavity 
functions in proximity of Rab becomes very small, they are subject to statistical errors. In 
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order to obtain smooth Duh-Haymet plots we needed to smooth the cavity functions obtained 
from the simulation. We did this by constructing the cubic smoothing spline which has as 
small a second derivative as possible. 



A. One component HS 

We carried out two simulations at pi ^ 0.650 (with a packing fraction of r/ = TcpiR\^/Q = 
0.340) and pi ^ 0.925 [r] = 0.484), the former corresponding to an intermediate density case 
and the latter to a liquid close to the freezing point. In our simulations we use Rn as unit 
of length. 

Inside the hard core, the bridge and the indirect correlation functions are monotonic and, 
for the cases here considered, there are no non-localities in the Duh-Haymet plots inside 
the core. Thus, to search for non-localities it is enough to analyze results in the external 
region. The resulting curves in the (-8,7) plane corresponding to points outside the hard 
core region are shown in Fig. ^ On the left the intermediate density case and on the right 
the high density one. We see that, as the density increases, the non-locality becomes more 
accentuated. Of course, the quality of a local approximation does depend on the choice 
of the correlation functions used as independent variable: plotting the bridge function as 
function of the direct correlation function we observed the opposite behavior. 

In order to compare the computer simulation results with the local approximate -8(7) 
relations used in the integral equations, we have plotted the curves corresponding to different 
closures: the hyper-netted chain (HNC) [l|: 

5(7) = , (12) 

the Percus Yevick (PY) Q; 

B{-f) = log{l + -f)-j , (13) 



the Martynov Sar 
Gazzillo (BPGG) 



dsov (MS) [13] and its generalization by Ballone, Pastore, Galli, and 

5(7) = (1 + "7)'/" -7-1, (14) 

(MS corresponds to a = 2, in the BPGG generalization a could be used as state dependent 
parameter to enforce thermodynamic consistence, here a fixed value of 15/8 has been used 



as suggested in and the modified Verlet (MV) jisl : 

BM = — - . (15) 

2(1 + 0.87] 

We can see that the best closures (MS, BPGG and MV), although not passing through the 
simulation curve, tend to follow its slope and curvature. When looking at Fig. ^one should 
also bear in mind that the values of the bridge function outside the hard core are not the 
most relevant for the quality of the structural and thermodynamic results of the closures. 



B. Equimolar binary mixture of AHS 

We carried out a simulation at pi = p2 — 0.589 [1] = 7r(pi_Rf^ + p2-R22)/6 = 0.375] and 
pi = 0.5. We chose -Rn = 1, Ru = 0.8, and R22 = 0.6. 

The results outside the hard core region are shown in the insets of the plots of Fig. |21 
There are non-localities in a neighborhood of the origin which corresponds to the large r 
region. These are more evident in the high density case. 

The most interesting feature shown in the figure is the difference between the curves at 
the two different densities. If the hypothesis of closures defined by a unique function -8(7) 
would be exact data for different densities should collapse into a unique curve in these plots. 
The data shown in Figs. ^ and |21 indicate clearly that this not strictly true. However, at 
low and intermediate densities the quantitative effect of the changing functional form is not 
dramatic. And even at the highest liquid densities, the success of closures such as MV, MS 
or BPGG can be probably explained in term of a higher sensitivity of the theory to localized 
(near the contact) features of the bridge functions more than to the behavior over the whole 
range of distances. 

Inside the hard core region the Duh-Haymet plots do not have non-localities. In Fig. El 
we show the results for the cavity functions yab for the system at the highest density. The 
plot for the unlike functions is more noisy than the plots for the like functions because yi2 
being smaller than yaa for a = 1, 2 is more subject to statistical error. 

In Fig. 121 we show the full Duh-Haymet plots for the system at the highest density, from 
the simulation (dots) and from integral equation theories (lines). The plots show how the 
MV approximation is the best one for this system. The unlike bridge function starts at 
r = close to the MV approximation, stays close to this approximation as r increases and 
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at some point have a smooth change in behavior and get closer to the PY curve. 

Fig. m shows the full bridge functions as a function of r for the system at the highest 
density. It is worth of notice the almost fiat region of the unlike bridge near the origin. 

C. Equimolar binary mixture of NAHS: Ru = R22, a < 

We carried out a simulation at pi = p2 — 0.573 (77 = 0.6). We chose Ru = R22 = 1 and 
R12 = 0.649 (a = —0.351). These radii values would be suitable for a reference system to 
model correlation in molten NaCl Q|. 

The results outside the hard core region are shown in the insets of the plots of Fig. |31 
There are non-localities in the neighborhood of the origin corresponding to the large r region. 

In Fig. iniwe show the results for the cavity functions yat- 

In Fig. we show the full Duh-Haymet plots from the Monte Carlo simulation (dots) 
and from the most common integral equation theories (lines). The approximation which 
seems to be closer to the like bridge function is MV: only at big r the bridge functions is 
well approximated by PY, MS, BPGG, and MV. The unlike bridge function starts at r = 
close to the PY approximation but as r increases it has a sudden change in behavior which 
displaces it away from all the approximations. Inside the hard core region the Duh-Haymet 
plots for the unlike functions exhibit significant non-localities in correspondence with the 
non monotonic behavior of the unlike cavity function (black dots in Fig. 

Fig. [7|shows the full bridge functions as a function of r. The unlike bridge function shows 
oscillations in a neighborhood of the origin. 

D. Equimolar binary mixture of NAHS: Ru = R22, a > 

We carried out a simulation at pi = p2 — 0.200 (77 = 0.209). We chose -Ru = R22 = 1 
and R12 = 1.2 {a = +0.2). Notice that this system undergoes phase separation when 
p = 2pi> 0.42. 

The results outside the hard core region are shown in the insets of the plots of Fig. |HI 
There are non-localities in a neighborhood of the origin corresponding to large distances. 

Also for this system, inside the hard core region the Duh-Haymet plots for the unlike 
functions have non-localities in a neighborhood of r = 0. These are smaller in extent than 
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the ones found for system C. In Fig. |Hlwe show the results for the cavity functions ijab- 

In Fig. ISl we show the full Duh-Haymet plots from the simulation (dots) and from 
the most common integral equations (lines). The approximations which seem to be closer 
to the like bridge function is MV and BPGG even if there is always a gap between the 
approximations and the simulation. The unlike bridge function starts at r = far away from 
all the approximations but as r increases it has a smooth change in behavior approaching 
the BPGG curve. 

Fig. ^1 shows the full bridge functions as a function of r. Again, the unlike bridge 
function have an almost flat behavior in a neighborhood of the origin. 

E. Equimolar binary mixture of NAHS: Rn ^ R22 

We carried out a simulation at pi = p2 — 0.589 {rj = 0.375). We chose -Rn = 1 and 
R12 = R22 = 0.6 (a = -0.2). 

The results outside the hard core region are shown in the insets of the plots of Fig. ITTl 
There are non-localities in a neighborhood of the origin which corresponds to the big r 
region. 

Inside the hard core region the Duh-Haymet plots have no non-localities. In Fig. IT2I we 
show the results for the cavity functions yab- 

In Fig. ^2 we show the full Duh-Haymet plots from the simulation (dots) and from the 
most common integral equations (lines). The approximation which is closer to the 11 bridge 
function is the MV. The one that is closer to the 22 bridge function is the BPGG. The 12 
bridge function starts at r = far away from all the 5 approximations and as r increases 
has a sudden change in behavior and starts following the BPGG approximation. 

Fig. El shows the full bridge functions as a function of r. The unlike bridge function 
shows again a qualitatively different behavior near the origin. 

IV. CONCLUSIONS 

From our analysis follows that the non-localities in the function relationship between the 
bridge functions and the indirect correlation functions may appear either outside of the hard 
core regions or inside of it. While the non-localities outside the hard core appear both in 
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the like and in the unhke functions, the ones inside the hard core appear only in the unlike 
functions (see Fig. Eland Fig. |HI), for the systems that we have studied. Their appearance 
can be directly related to the peculiar behavior of the unlike cavity correlation function 
inside the hard core. 

As is shown by a comparison of the plots of Fig. ^ and from Fig. El the non-localities 
become more accentuated as we increase the coupling (the density) of the system. Nonethe- 
less Fig. ISl shows that the non-localities may appear even in a weakly coupled system (in 
this case symmetric NAHS with positive non additivity). Among the systems studied the 
one which presents the worst non-localities is the equimolar symmetric NAHS with negative 
non additivity (see Fig. . For this system the Duh-Haymet plot for the unlike bridge 
function is non-local both in the hard core region (in a neighborhood of r = 0) and outside 
of it (at large r). 

We can conclude that the two hypothesis of a local function approximation for the bridge 
functionals of the indirect correlation functions and the stronger hypothesis of unique func- 
tional form independent on the state, are not strictly supported by the numerical data. For 
the one component system, this finding is consistent with the observed density dependence 
of the bridge function reported in [17]. We observe clear violations of both the assumptions 
increasing with the density. This negative statement should be somewhat mitigated by re- 
alizing that the violations of the locality, in the systems studied, are limited to the small 
and large distances regions. The latter, corresponding to the region of the fast vanishing 
of the bridge functions affect very little the thermodynamic and structural properties of 
the systems. The former are presumably more important for the level of thermodynamic 
consistence of the theory but have small effect on quality of the structural results. The well 
known success of closures like MS, BPGG and MV supports such point of view. 

From comparison with the simulation data in the cases we have studied, we conclude 
that the best approximations of the true hard sphere bridge functions are provided by the 
MV and BPGG even if, especially in the unlike bridge functions, there are a wide variety of 
characteristic behaviors which are not captured by any of the most popular integral equa- 
tion approximations. In this respect, we feel that a final comment on the local functional 
approximation in the case of multicomponent systems is in order. Indeed, density functional 
theory allows to say that the bridge function Bij should be a functional of all the pair cor- 
relation functions, not only the one. Thus, we could have a function approximation 



12 



-6ij(7ii(^)) 7i2('"), 722('")) which would be local in space but not with respect to the com- 
ponents. At the best of our knowledge, up to now no attempt has been done to explore 
this additional freedom to improve the modeling of the bridge functions in multicomponent 

systems. 
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LIST OF FIGURES 

Fig. The first two graphs are Duh-Haymet plots (dots), outside the hard core region, for the 
one component HS system (the lines show the behavior of integral equation closures). 
On the left p = 0.650 on the right p = 0.925. 

Fig. 121 Full Duh-Haymet plots obtained by the inversion of the Monte Carlo simulation data 
(dots) compared with some of the most common integral equation theories (lines) for 
the equimolar binary mixture of AHS at two different densities (in the second and 
third plot only results at the highest density are shown). Ru = 1, = 0.8, and 
R22 = 0.6. The insets shows the portion of the bridge function outside the hard cores. 

Fig. 121 Cavity functions inside the hard core for the equimolar binary mixture of AHS (at the 
same conditions as in Fig. |21at the highest density). The plot shows the behavior of 
the functions defined in © and ^ (notice the logarithmic scale on the ordinates), 
the triangles denote the 22 function, the open circles the 11 function, and the closed 
circle the 12 function. 

Fig. m Bridge functions Bab{r) for the equimolar binary mixture of AHS (at the same con- 
ditions as in Fig. |21 at the highest density). The insets shows magnifications of the 
regions just outside of the hard cores. 

Fig. El Full Duh-Haymet plots obtained by the inversion of the Monte Carlo simulation data 
(dots) and by some of the most common integral equation theories (lines) for the 
equimolar binary mixture of NAHS with equal like diameters and negative non addi- 
tivity a = —0.351, at pi = 0.589. Ru = R22 = 1 and R12 = 0.649. The insets shows 
the portion of the bridge function outside the hard cores. 

Fig. ini Cavity functions for the equimolar binary mixture of NAHS with equal like diameters 
and negative non additivity (at the same conditions as in Fig. E}. The graph shows 
the behavior of the functions defined in © and (notice the logarithmic scale on the 
ordinates), the open circle denotes the like functions and the closed circle the unlike 
one. 

Fig. [71 Bridge functions Bab{r) for the equimolar binary mixture of NAHS with equal like 
diameters and negative non additivity (at the same conditions as in Fig. EJ. The 
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insets shows magnifications of the regions just outside of the hard cores. 

Fig. IHl Full Duh-Haymet plots obtained by the inversion of the Monte Carlo simulation data 
(dots) and by some of the most common integral equation theories (lines) for the 
equimolar binary mixture of NAHS with equal like diameters and positive non addi- 
tivity a = +0.2, at pi = 0.200. -Rn = R22 = 1 and R12 = 1.2. The insets shows the 
portion of the bridge function outside the hard cores. 

Fig. El Cavity functions for the equimolar binary mixture of NAHS with equal like diameters 
and positive non additivity (at the same conditions as in Fig. The graph shows the 
behavior of the functions defined in © and (jHI) (notice the logarithmic scale on the 
ordinates), the open circle denotes the like functions and the closed circle the unlike 
one. 

Fig. El Bridge functions Bab{r) for the equimolar binary mixture of NAHS with equal like 
diameters and positive non additivity (at the same conditions as in Fig. The 
insets shows magnifications of the regions just outside of the hard cores. 

Fig. ^2 Full Duh-Haymet plots obtained by the inversion of the Monte Carlo simulation data 
(dots) and by some of the most common integral equation theories (lines) for the 
equimolar binary mixture of NAHS with different like diameters Ru = 1 and R12 = 
R22 = 0.6, at pi = 0.589. The insets shows the portion of the bridge function outside 
the hard cores. 

Fig. ^1 Cavity functions for the equimolar binary mixture of NAHS with different like diam- 
eters (at the same conditions as in Fig. [TT|) . The graph shows the behavior of the 
functions defined in (jH)) and © (notice the logarithmic scale on the ordinates), the 
triangles denote the 22 function, the open circles the 11 function, and the closed circle 
the 12 function. 

Fig. ^1 Bridge functions Bab{r) for the equimolar binary mixture of NAHS with different like 
diameters (at the same conditions as in Fig. fTT|l . The insets shows magnifications of 
the regions just outside of the hard cores. 
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FIG. 2: R. Fantoni and G. Pastore 
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FIG. 3: R. Fantoni and G. Pastore 
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FIG. 8: R. Fantoni and G. Pastore 
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FIG. 10: R. Fantoni and G. Pastore 
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